Collective dynamical response of coupled oscillators with any network structure 
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We formulate a reduction theory that describes the response of an oscillator network as a whole 
to external forcing applied nonuniformly to its constituent oscillators. The phase description of 
multiple oscillator networks coupled weakly is also developed. General formulae for the collective 
phase sensitivity and the effective phase coupling between the oscillator networks are found. Our 
theory is applicable to a wide variety of oscillator networks undergoing frequency synchronization. 
Any network structure can systematically be treated. A few examples are given to illustrate our 
theory. 
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An assembly of coupled limit-cycle oscillators often be- 
haves like a single large oscillator. This general scenario 
recurs in a wide variety of rhythmic phenomena in living 
organisms, ranging from circadian oscillations, cardiac 
rhythms to pathological phenomena such as epilepsy and 
Parkinsonian disease [J, H, 0, 0, El . Recent experiments 
using electrochemical oscillators simulate such naturally 
arising populations of oscillators in an idealized form [6| . 

Many previous studies have been devoted to answer- 
ing how and under what conditions oscillators mutually 
synchronize. In comparison, little attention has been 
paid to investigating the dynamical response of an oscil- 
lator network to external stimuli. Such an inquiry would 
shed light on mechanisms underlying biological functions 
(such as the phase response curves of circadian rhythms 
7]), external control, and inter-network synchronization 
of oscillator networks [8|. Establishing a description of 
the collective dynamics on the basis of "microscopic" 
knowledge, i.e., the nature of the constituent oscillators 
and their mutual coupling, presents a challenging task 
from a theoretical point of view. In particular, it would 
be ideal if the collective dynamics could be described in 
terms of a single, suitably defined collective phase in a 
closed form. A quantity of central importance will then 
be the phase sensitivity of a population as a whole to ex- 
ternal weak stimuli 14]. This function has already been 
derived [9( for a system consisting of a large assembly of 
identical phase oscillators with global coupling, with each 
oscillator being driven independently by random time- 
dependent noise. 

In this Letter, we argue that there is yet another gen- 
eral situation for which a similar phase description can 
be formulated. In contrast to Ref. 0, we consider noise- 
free but nonidentical oscillators undergoing full frequency 
synchronization (in which all the oscillators have exactly 
the same frequency due to coupling) . A similar situation 



has also been studied by Ko and Ermentrout for two 
symmetrically coupled and globally coupled oscillators 
very recently [lj| • A distinct advantage of our present 
approach is that it may deal with any system size, any 
connectivity, any heterogeneous coupling, and nonuni- 
form external forcing. Moreover, the theory is extended 
to include multiple populations simply by reinterpreting 
the external stimuli applied to a given population as the 
coupling forces originating from the other populations, 
which enables us to predict the synchronization behav- 
ior between oscillator networks. General formulae for the 
collective phase sensitivity and the effective phase cou- 
pling between the oscillator networks are found. 

Consider a network of N coupled limit-cycle oscillators 
under external forcing. The forcing is generally nonuni- 
form, i.e., individual oscillators receive different inputs. 
As is well known Q, if the heterogeneity of oscillators, 
the coupling between oscillators, and external forcing are 
weak, the system is describable by the phase equation 



N 



= uji + ^T ij ((j) i -4>j) +eZ(d>i)£i(t). 



(1) 



Here fa is the phase of the ith oscillator (i = 1, . . . , N), 
LOi its natural frequency, and IY7 the coupling force from 
the jth oscillator to the ith oscillator. The terms 
and Z(fa) respectively represent the time-dependent ex- 
ternal force and the phase sensitivity of the oscillator i 
to external perturbation [15[ . Parameter e is the charac- 
teristic intensity of the external forcing. 

Our aim is to establish the collective phase description 
for Eq. ([T]), i.e. to derive the dynamical equation for a 
suitably defined macroscopic variable that describes the 
response to external forcing. This is generally formulated 
under two basic assumptions, (i) In the absence of ex- 
ternal forcing a stable periodic solution corresponding to 
a fully frequency-synchronized state exists, and thus, the 
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oscillator network behaves as a single large limit-cycle 
oscillator, (ii) The external force is even weaker than 
the coupling force, i.e. e <C 1, so that the synchronized 
state is almost unaltered under external forcing. Under 
these assumptions, the phase reduction method Q ap- 
plicable to a weakly perturbed oscillator can be applied 
(once again) to the oscillator network by interpreting the 
unperturbed system as a single limit-cycle oscillator. 

For convenience, we begin by rewriting Eq. (TTJ) in terms 
of the A-dimensional state vector X = [<j>\, <fo, ■ ■ ■ , 4>n) 
as 



X=F(X) + ep(X,t), 



(2) 



where Fi(X) = lo { + ~ 4>j) and Pi{ x ,t) = 

Z (<f>i)£i(t) . A frequency synchronized state (or, more ex- 
actly, a phase locked state) for e = is found as a solution 
of Fi(X) = fl for all i. This "limit-cycle" solution is de- 
noted by 

Xo(e) = (e + 0°,e + </>",..., e + ^V), e = nt, (3) 

where </>° is constant. We then extend the definition of 
outside the limit-cycle orbit as a scalar field Q(X). The 
identity 6 = (dQ/dX) ■ {dX/dt) implies 



e = 



do 

dX 



{F(X) + ep(X,t)}. 



(4) 



As is usually done in the phase reduction process, we now 
adopt the definition of Q(X) such that in the absence of 
the forcing, 8 = f2 is identically satisfied, i.e., (dQ/dX) ■ 
F(X) = Q. We call Q(X) the collective phase Due 
to the assumption of weak forcing, (dO/dX) ■ p(X,t) 
may be evaluated at X = Xq(Q). Then, to the lowest 
order in e Eq. (j4]) becomes 



e = n + e^-- P (x a ,t). 



(5) 



Let the linearized equation of ([2]) with e = be Y = LY, 
where Y is the deviation defined by Y = X — X . Due 
to the symmetry of F(X), the Jacobian L is a constant 
matrix, one eigenvalue of which equals zero with the cor- 
responding eigenvector given by U = dX (Q)/dQ = 
(1, 1, . .., 1). We define row vector U* as the left zero- 
eigenvector of L, i.e., U*L = 0, with the normalization 
condition U*U=1, or, J^Zi u * = L :t 

can then be ar- 
gued that dQ/dXo becomes identical to U* [2]. Thus, 
Eq. ([5]) takes the form 



e = n + c c(e) •€(*), 



(6) 



where £ = (£i, £2, ■ ■ • , £,n) and C(0), which is interpreted 
as the collective phase sensitivity, is a vector with com- 
ponents 



Ci{G) = UTZ(Q + </,° i ). 



(7) 



Equations (O and ([7]) have clarified the response of the 
collective mode of the synchronized network. There, the 
forcing to the oscillator i has turned out to be weighted 
by U* . We thus call U* the weight vector henceforth. 

In what follows, we often consider uniform external 
forcing, £i(t) = for all i. In such a case, Eq. ^ is 
further reduced to = ft + e£(0)£(i), where C(0) is the 
collective phase sensitivity defined for uniform forcing, 
C(©) = E,=i U*Z(S + Generally speaking, ({&) 

deviates from Z(<j>) more significantly as the phases of 
the constituent oscillators are more widely distributed. 

An analytic formula of weight vector U* is given as 
follows. The elements of L are given by 



N 



(8) 



Note that the relation LU = L%j — holds. 

We define the (i, i)-cofactor of L as Mj = det L(i,i), 
which is the determinant of the submatrix L(i,i), that 
is L with the i-th row and column removed. One can 



prove that, for any L with Ylj=i -^ij — 0> the row vec- 
tor (Mi, M2, Mn) is a left zero-eigenvector of L, i.e., 
Y^jLi MjLji = [IH. Using the normalization condi- 
tion, we obtain the algebraic expression of U* as 



U* = (M 1 ,M 2 ,...,M N )/M, 



(9) 



where M = Mj. This expression is valid for any 

network structure. Note that, for networks with both 
bidirectional connection patterns and symmetric cou- 
pling functions, L is symmetric and U* is trivially l/N 
for all i, which is the case even for networks with strongly 
heterogeneous connectivity including scale-free networks. 
In many cases, however, L is asymmetric and U* is het- 
erogeneous. 

Next, we formulate the collective phase description of 
multiple networks of phase oscillators. We are concerned 
with the case in which external forcing is absent, while 
e£j(t) in the last term in Eq. ([T]) is interpreted as rep- 
resenting the coupling force coming from oscillators of 
other networks. For clarity, we consider a simple system 
in which two identical networks composed of N oscilla- 
tors, called group A and group B, are uniformly cou- 
pled (the extension to a more general case, e.g., weakly 
heterogeneous multiple groups, is straightforward). The 
dynamical equations of the system are given by 



JV 



# = Ui + Tiiitf ~ <i>f) + eZ^mk}), 



N 



(10) 



tf) + eZ(tf )£({<%}) 



where <ff is the phase of the oscillator i in group X (X=A, 
B), and £({</>£}) denotes a function of (j>f, . . . , 4>n, 
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which represents the uniform coupling force coming from 
group X. Denoting the collective phases of the respective 
groups by 9a and &b, we obtain the resulting phase 
equation in the form 

e A = n + eC(e A )£({e B + ^}), 

e B = tt + eC(e B )£({e A + ^}). ' 

Noting that the dominant time-dependence of 0a, b is 
fit, we may obtain the effective coupling 7(0a — ©b) 
between the groups by time-averaging of Eq. (fTTj) over 
the common period 2ir/Q. Putting Oa.b = ^ + $a,b, we 
perform the averaging as 

7(e A - e B ) = 7(#a - B ) = 
-L f * d(nt)c(n* + (>A)tmt + e B + (12) 

2tt Jo 

where 7 rather than T has been used to indicate that this 
coupling function acts between the groups. In this way, 
we succeeded in deriving the collective phase equation in 
the simple form 

e A = tt + e7(e A -e B ), 

e B = n + e 1 (Q B ~Q A ). 

We give two examples to illustrate our theory. 

Example I: the dependence of weight vector U* on 
the connectivity. We consider a network of identical 
oscillators with homogeneous coupling, i.e., u>i = lu 
and Tij(A(f>) = AyT(A0) undergoing perfect phase 
synchrony. Here, A is the adjacency matrix describ- 
ing the connectivity, which is generally asymmetric and 
weighted. In such a system, the weight vector depends 
solely on the network architecture. By rescaling t and e, 
we put r'(0) = — 1 without loss of generality. In perfect 
phase synchrony, i.e., = (fP for all i, L is a Lapla- 
cian matrix generalized for asymmetric and weighted net- 
works, given by = -Aij for i ^ j and Lu = Ysj^a Ar 
We consider two small networks in which the weight vec- 
tor U* is easily calculated via the algebraic expression, 
Eq. ((9]). Figure QJa) is a weighted network, where the 
weight vector is found to be a simple reflection of the 
connection weights. Figure [ljb) is a non-weighted net- 
work but its adjacency matrix is asymmetric. It is worth 
noticing that oscillator 2 is more influential than oscil- 
lator 3, although they have locally the same topological 
properties: one inward and two outward connections. In 
general, the weight vector depends on the global topol- 
ogy 

Example II: collective phase sensitivity and group syn- 
chronization in limit-cycle oscillators. We illustrate that 
the collective phase sensitivity of a group of coupled oscil- 
lators varies with intra-group coupling strength. We then 
consider two groups of coupled oscillators with an addi- 
tional inter-group coupling of fixed strength, and show 




FIG. 1: Examples of weight vector U*. Using Eq. (J5J), we find 
(a) U{ :UZ = 1:k and (b) U{ : f/ 2 * : [/■? : UX = 2 : 4 : 3 : 1. 

that a nontrivial qualitative change in the synchroniza- 
tion behavior between the groups occurs when the indi- 
vidual collective phase sensitivities changes as a result of 
modifications to the intra-group coupling strength. 

As schematically illustrated in Fig. [DJa), we consider 
the system in which a pair of identical groups A and B, 
each of which consists of two coupled limit-cycle oscilla- 
tors, are mutually coupled. We use the Hindmarsh-Rose 
model as the limit-cycle oscillator, a model originally pro- 
posed as a neural model. The system reads 

4 

±i = Zxt-Xj+yj-fii+y^^DijXj, yi = \-hx\-yi. (14) 

Here, the coupling matrix D is given as C12.21.34.43 = 
K, D n , 22,33,44 = ~K,Di3, 23,31,41 = e with K and e be- 
ing the coupling intensities for intra- and inter-groups, 
respectively. We assume 1 > X > e = 1.0 x 10 -5 . We 
set £*i = /i3 = —3.000 and [12 = /M = —3.001, corre- 
sponding to 0^1,3 ~ 1.80476,072,4 — 1.80443, and thus 
Auj ~ 3.3 x 10~ 4 . The wave form x((p) and the phase 
sensitivity Z(ifi) of an isolated oscillator obtained numer- 
ically are shown in Figs. [2jb) and [2c), respectively. 

The corresponding phase model is given by Eq. (fTOf . 
where T y (A0) = KT(A^), £({0£}) = xfa), and 
= x {4>?>)- From the phase reduction theory, 
the coupling function is calculated as r(</>i — </> 2 ) = 
i f 2 * d(u)t)Z{<j> x + u)t){x(4> 2 + uit) - xfa + ujt)}. For 
convenience, we display r a (A0) = T(A<j)) — T(— A<fi) in 
Fig. HJd). The phase difference Acffi between the oscilla- 
tors 1 and 2 of a synchronized state is found as a stable 
solution of 4>\ — 4>2 (where e = is assumed), and thus, a 
solution of r a (A</>°) = Auj/K. The predicted phase dif- 
ference is plotted in Fig. [^e) as a curve. It agrees well 
with numerical data obtained through direct numerical 
integration of Eq. flU]). Using Z((p), AcjP , T(A<j)), and its 
derivative T'(A(f>) obtained numerically, we can calculate 
Ui,U%, and C(©)- The results are shown in Fig. ^i) and 
its caption. For large K (compared to Auj), C(©) is indis- 
tinguishable from Z ((p). As K decreases, AcjP becomes 
larger, resulting in considerable variation in C(@)- 

Given C(@)> the synchronization behavior between 
groups is now predicted. The collective coupling func- 
tion 7(A8) is calculated from Eq. (fT2|) where = x(<f>) 
in the system under consideration. For convenience, we 
display the antisymmetric part of 7(A0) in Fig. [2tg). 
Putting G>a = <9b, we find the stable phase locking so- 
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broad applicability would be expected. 
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FIG. 2: (color online) Results for a network of limit-cycle 
oscillators, (a) Network structure under consideration, (b) 
Waveform x((j>) and (c) phase sensitivity Z(<f>) of an indi- 
vidual oscillator, (d) Antisymmetric part of the coupling 
function, r a (A0). (e) Phase difference A</> between oscil- 
lators in each group, (f) Collective phase sensitivity C(9) 
defined for uniform forcing. For K = 0.007, 0.002, 0.0015, 
(Ui,US) is respectively about (1.35, -0.35), (2.87, -1.87) and 
(2.89, —1.89). (g) Antisymmetric part of the collective cou- 
pling function, 7 a (AO). (h) Phase difference AO between the 
groups. 



lution between groups, AO = 0a — 0b- Predicted A0 
is exhibited by the curve in Fig. H^h), implying that the 
in-phase solution becomes unstable at around K — 0.006 
(via a pitch-fork bifurcation) and the out-of-phase so- 
lution appears below. Phase difference A0 (or equiva- 
lently, <f>i— 03 ) obtained from direct numerical integration 
of Eq. ([HI) is plotted in Fig. H^h), which convinces us of 
the precision of the phase description. 

In summary, we have formulated the response of the 
collective phase to weak external forcing given to con- 
stituent oscillators and the phase description for inter- 
acting oscillator networks. The present theory is valid 
for a wide class of weakly coupled oscillator networks 
undergoing full frequency synchronization, and thus, a 
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The phase sensitivity is the function of the phase of an 
oscillator describing the phase shifts (compared to unper- 
turbed one) normalized by the intensity of weak pulsative 
perturbation given to the oscillator 0, @, 0]. It is also 
sometimes called the infinitesimal phase response curve 

m- 

When the forcing is given by a vector £j(t) in the orig- 
inal dynamical system model before reduction, the cor- 
responding forcing term in the phase-reduced equation 
is generally given, unlike Eq.{T}, by a scalar product be- 
tween this forcing vector and the phase sensitivity vector 
Z(4>i) of the ith oscillator Q]. Here, Z(cj>i) is the left 
eigenvector with the zero eigenvalue defined for the lin- 
earized system of the ith oscillator about its limit-cycle 
orbit. However, to avoid unnecessary complications, we 
have used a simpler form of the forcing term in Eq. Jl]) 
assuming that the original forcing vector has only one 
nonvanishing component fj. 

In our particular system, O can be explicitly given for 
X ~ Xo as O = U* X + ©o with some constant 9o, 
because 9 = U*X = U*F(X )+U*LY = J^Zi U i n = 
H. Correspondingly, one may immediately obtain Eq. (|6]> 
by applying U* to the lowest order description of Eq. ([2j, 
X = F(X ) + LY + ep(X ,t), as well as by following 
the normal phase reduction process. 



